
/*!
 * @file normal.h
 * @ingroup MathUtil
 *
 * Routines for generating random numbers.
 *
 *
 */



#ifndef _QT_NORMAL_H_
#define _QT_NORMAL_H_

#include <cstdlib>
#include <cmath>
#include <cassert>

#ifndef M_TWOPI
#define M_TWOPI 6.2831853071795864769L
#endif

#ifndef M_SQRTPI
#define M_SQRTPI 1.7724538509055160273L
#endif


#include <boost/random.hpp>

#include "mathutil.h"
#include "randomutil.h"

namespace qt {

namespace random {

/*!
 * This class provides random deviates from a
 * normal (0,1) distribution using the polar transform method.
 */
template <typename UniformGenerator >
class PolarStdNormalGenerator {

public:
    typedef double result_type;


private:

    PolarStdNormalGenerator();
    //      PolarStdNormalGenerator(const PolarStdNormalGenerator &);


public:

    PolarStdNormalGenerator(UniformGenerator & gen)
        : _cached(false), _ugen(gen) {}

    double operator()() {
        if(_cached) {
            _cached = false;
            return _cache;
        }
        double x1 = _ugen();
        double x2 = _ugen();
        double y1 = std::sqrt(-2*std::log(x1))*std::cos(M_TWOPI*x2);
        double y2 = std::sqrt(-2*std::log(x1))*std::sin(M_TWOPI*x2);
        _cache = y2;
        _cached = true;
        return y1;
    }

    void seed(unsigned long s) {
        _ugen.seed(s);
        _cached = false;
    }

private:
    bool _cached;
    double _cache;
    UniformGenerator & _ugen;
};



/*!
 * This class provides access to standard functions
 * associated with the standard normal (Gaussian) distribution.
*/
class StdNormal {

public:
    /*!
     * Cumulative normal distribution.
     *
     * @param x is a real number
     * @returns Prob(X<x) where X is distributied N(0,1).
     *
     * History:
     *
     */
    static double cdf(double x) {
        return 0.5*qt::mathutil::erfcc((-1.0/M_SQRT2)*x);
    }

    /*!
     * Inverse of the cumulative normal distribution function.
     *
     * @param u a real number strictly between 0.0 and 1.0
     * @returns x the largest number such that P(X<x) < u.
     *
     * History:
     *
     */
    static double invCdf(double u) {
        return M_SQRT2*qt::mathutil::erfinv(2*u-1);
    }


    /*!
     * The density function of the Normal distribution with respect to dx.
     *
     * @param x a real number
     * @returns Prob( x-dx < X < x+dx )/dx
     *
     * History:
     *
     */
    static double pdf(double x) {
        return (1.0/(M_SQRT2*M_SQRTPI))*std::exp(-0.5*x*x);
    }



    template <typename GeneratorType>
    class GeneratorWrapper {

    public:

        typedef PolarStdNormalGenerator<GeneratorType> Generator;

        typedef WrappedGenerator
        < GeneratorType, PolarStdNormalGenerator<GeneratorType> >
        BoundGenerator;

    };


    typedef
    GeneratorWrapper<DefaultGeneratorType>::Generator DefaultGenerator;

    typedef
    GeneratorWrapper<DefaultGeneratorType>::BoundGenerator
    DefaultBoundGenerator;


    /*!
     * A factory method.
     */
    static DefaultGenerator getDefaultGenerator() {
        qt::random::DefaultGeneratorType & dgen =
            qt::random::DefaultGeneratorType::getDefaultGenerator();

        DefaultGenerator gen(dgen);
        return gen;

    }

    template <typename GeneratorType>
    static typename GeneratorWrapper<GeneratorType>::Generator
    getGenerator(GeneratorType & ugen) {
        typename GeneratorWrapper<GeneratorType>::Generator gen(ugen);
        return gen;
    }

    /*! Returns an instance of a generator that owns
     *  its own base generator.  Reseeding this generator
     *  will not affect other generators.
    */
    template<typename GeneratorType>
    static
    typename GeneratorWrapper<GeneratorType>::BoundGenerator
    getBoundGenerator(unsigned long seed = default_seed) {

        typename GeneratorWrapper<GeneratorType>::BoundGenerator wgen;

        wgen.seed(seed);

        return wgen;
    }

};

/*!
 * This is a workaround for some issues with the compiler not instantiating
 * templates in the right spot.
 */
class NormalBase {
public:
    /*!
     * Cumulative distribution.
     *
     * @param x is a real number
     * @returns Prob(X<x) where X is distributied U(0,1).
     *
     * History:
     *
     */
    static double cdf(double x, double mu, double sigma) {
        return StdNormal::cdf((x-mu)/sigma);
    }

    /*!
     * Inverse of the cumulative distribution function.
     *
     * @param u a real number strictly between 0.0 and 1.0
     * @returns x the largest number such that P(X<x) < u.
     *
     * History:
     *
     */
    static double invCdf(double u,double mu, double sigma) {
        assert(u>=0 && u<=1);
        return StdNormal::invCdf(u)*sigma + mu;
    }


    /*!
     * The density function of the  distribution with respect to dx.
     *
     * @param x a real number
     * @returns Prob( x-dx < X < x+dx )/dx
     *
     * History:
     *
     */
    static double pdf(double x,double mu, double sigma) {
        return StdNormal::pdf((x-mu)/sigma);
    }

private:
    double _mu;
    double _sigma;

public:
    /*! default constructor */
    NormalBase(double mu=0.0, double sigma=1.0) : _mu(mu), _sigma(sigma) {}

    double get_mu() {
        return _mu;
    }

    double get_sigma() {
        return _sigma;
    }

    void set_mu(double mu) {
        _mu = mu;
    }

    void set_sigma(double sigma) {
        _sigma = sigma;
    }

    /*!
     * Cumulative distribution.
     *
     * @param x is a real number
     * @returns Prob(X<x) where X is distributied N(mu,sigma).
     *
     * History:
     *
     */
    double cdf(double x) {
        return NormalBase::cdf(x,_mu,_sigma);
    }

    /*!
     * Inverse of the cumulative distribution function.
     *
     * @param u a real number strictly between 0.0 and 1.0
     * @returns x the largest number such that P(X<x) < u.
     *
     * History:
     *
     */
    double invCdf(double u) {
        return NormalBase::invCdf(u,_mu,_sigma);
    }


    /*!
     * The density function of the  distribution with respect to dx.
     *
     * @param x a real number
     * @returns Prob( x-dx < X < x+dx )/dx
     *
     * History:
     *
     */
    double pdf(double x) {
        return NormalBase::pdf(x,_mu,_sigma);
    }


};

/*!
 * This class provides access to standard functions
 * associated with the normal(mu,sigma) distribution, where
 * mu is the mean and sigma is the standard deviation.
 */
class Normal : public NormalBase {


public:
    /*! default constructor */
    Normal(double mu=0.0, double sigma=1.0) : NormalBase(mu,sigma) {}


    /* a helper nested class */
    template <typename GeneratorType>
    class GeneratorAggregateRef {
    private:
        NormalBase _normal;
        typename StdNormal::GeneratorWrapper<GeneratorType>::Generator _ngen;

    public:
        GeneratorAggregateRef(double mu, double sigma,
                              GeneratorType & gen)
            : _normal(mu,sigma), _ngen(gen) {}

        double operator()() {
            return _normal.get_mu() + _normal.get_sigma()*_ngen();
        }

    };

    /* a helper nested class */
    template <typename GeneratorType>
    class GeneratorAggregateBound {

    public:
        GeneratorAggregateBound(double mu, double sigma)
            : _normal(mu,sigma), _ugen(), _ngen(_ugen) {}


        double operator()() {
            return _normal.get_mu() + _normal.get_sigma()*_ngen();
        }

        void seed(unsigned long x=default_seed) {
            _ugen.seed(x);
        }

    private:

        GeneratorType _ugen;
        NormalBase _normal;

        typename
        StdNormal::GeneratorWrapper<GeneratorType>::Generator _ngen;

    };





    template <typename GeneratorType>
    class GeneratorWrapper {

    public:

        typedef GeneratorAggregateRef<GeneratorType> Generator;
        typedef GeneratorAggregateBound<GeneratorType> BoundGenerator;


    };

    typedef
    GeneratorWrapper<DefaultGeneratorType>::Generator DefaultGenerator;
    typedef
    GeneratorWrapper<DefaultGeneratorType>::BoundGenerator DefaultBoundGenerator;

//       typedef
//       GeneratorWrapper<DefaultGeneratorType>::BoundGenerator
//       DefaultBoundGenerator;

    /*!
      * A factory method.
      */
    static
    DefaultGenerator
    getDefaultGenerator(double mu, double sigma) {

        DefaultGeneratorType & dgen =
            DefaultGeneratorType::getDefaultGenerator();


        DefaultGenerator gen(mu,sigma,dgen);
        return gen;

    }

    template <typename GeneratorType>
    static typename GeneratorWrapper<GeneratorType>::Generator
    getGenerator(GeneratorType & ugen, double mu, double sigma) {

        typename
        GeneratorWrapper<GeneratorType>::Generator gen(mu,sigma,ugen);

        return gen;
    }


    /*! Returns an instance of a generator that owns
     *  its own base generator.  Reseeding this generator
     *  will not affect other generators.
    */
    template<typename GeneratorType>
    static
    typename GeneratorWrapper<GeneratorType>::BoundGenerator
    getBoundGenerator(double mu, double sigma,
                      unsigned long x = default_seed) {

        typename GeneratorWrapper<GeneratorType>::BoundGenerator
        gen(mu,sigma);
        gen.seed(x);
        return gen;

    }



};

}; // namespace qt::random

}; // namespace qt::


#endif
